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Oh Abstract 
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An algorithm for perfect simulation from the unique solution of the 

distributional fixed point equation Y —d UY + U(1 — U) is constructed, 

^ where Y and U are independent and U is uniformly distributed on [0, 1] . 

This distribution comes up as a limit distribution in the probabilistic 

f^i analysis of the Quickselect algorithm. Our simulation algorithm is 

^. based on coupling from the past with a multigamma coupler. It has 

^N four lines of code. 
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1 Introduction 



In a probabilistic analysis of the algorithm Quickselect Hwang and Tsai [8] 
^ showed that, when applied to a uniformly random permutation of length n 

and selecting a rank of order o(n), the normalized number of key exchanges 
performed by Quickselect converges in distribution to a limit distribution 
(j,. This limit distribution is characterized as the unique probability measure 
fj, = C(Y) such that 

Y = UY + U(1-U), (1) 
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where = (also =d) denotes equality in distribution and U is uniformly dis- 
tributed over the unit interval [0, 1] and independent of Y. 

The distribution \i was studied in [9]. In particular we showed that /x 
has a bounded, l/2-H61der continuous density, (j, is supported by the unit 
interval [0, 1] and we developed a method to numerically approximate the 
density and the corresponding distribution function. In Remark 2.9 of [9] 
we noted that this is sufficient to theoretically construct an algorithm for 
perfect simulation from /x based on von Neumann's rejection method along 
the approach taken in Devroye [2]. While the numerical approximations 
yield an algorithm for perfect simulation in almost surely finite time, the 
convergence rates of our approximations are poor and the expected running 
time is infinite. We do not expect such an algorithm to terminate within 
our lifetimes. 

Recently, Fill and Huber [6j published an algorithm for perfect simu- 
lation of a related distribution, known as the Dickman distribution and 
characterized as unique solution of the distributional fixed point equation 
Y =d UY+1. This algorithm is based on coupling from the past of a Markov 
chain with the Dickman distribution as stationary distribution. The method 
makes use of a multigamma coupler and of a dominating chain to deal with 
the unbounded support of the Dickman distribution. In fact Fill and Huber 
develop their algorithm for a more general class of distributions, the Vervaat 
perpetuities. Devroye and Fawzi [3] presented a different multigamma cou- 
pler and a different dominating chain resulting in a faster coupling from the 
past algorithm for the Dickman distribution. Both algorithms are also fully 
satisfactory from a practical point of view, millions of independent samples 
from the Dickman distribution can be generated within seconds. 

In this note we construct a coupling from the past algorithm for the 
solution fi of M. Compared to the more difficult Dickman case we benefit 
from the special analytic structure of the densities <p x of Ux + U(1 — U) for 
x £ [0, 1]. In particular, we have 

inf inf <p x (t) > 1/2, (2) 

ze[o,i]te[o,i/4] 

which allows for the construction of a multigamma coupler as proposed by 
Murdoch and Green |l(Jl Section 2.1]. This results in a fast and simple 
four-line-code algorithm. 

Note that a general method described in an unpublished extension of |3] , 
see Fawzi [5], can also be applied to our fj,: In [5, Section 4] it is shown that 
when one is able to perfectly simulate from the solution of Y =d AY + 1 
with a random < A < 1 this can be turned into an algorithm to simulate 



from the solution of Y =d AY + B, whenever B > is bounded. Here, 
(A, B) is independent of Y. Hence, this method together with the simula- 
tion algorithm for the Dickman distribution yields as well an algorithm to 
simulate from /x. 

For general perfect simulation algorithms for another class of perpetuities 
see Devroye and James |3]. For perfect simulation algorithms from station- 
ary distributions of positive Harris recurrent Markov chains see Hobert and 
Robert [7J. 

In the field of exact simulation from nonuniform distributions it is cus- 
tomary to assume that a sequence of independent and identically, uniformly 
on [0, 1] distributed random variables is available and that elementary op- 
erations of and between real number such as +, — , /, *, yfx, logx, etc., can 
be performed with absolute precision, see Devroye pQ for a comprehensive 
account on nonuniform random number generation. 

2 Markov chain and multigamma coupler 

An underlying ergodic Markov chain (Xj) on [0, 1] having fj, as stationary 
distribution is given as follows: For all x £ [0, 1], given Xj = x, we define 
-Xj.fi to be distributed as Ux + U(l — U) with a uniform [0, 1] random 
variable U. In the context of coupling from the past a realization of such 
a Markov chain is usually constructed with a deterministic update function 
$ : [0, 1] x [0, 1] -)■ [0, 1] such that X j+1 := $(Xj, U j+ i) yields a realization 
of the chain, where (Uj) is a sequence of independent and uniform [0, 1] 
random variables. A trivial choice for <& is (x, u) h- >■ ux + u(l — u). However, 
to make coupling of the chains possible, we follow the construction of a 
multigamma coupler as described by Murdoch and Green |10j . 

The construction is as follows: Assume that a probability density / 
is written as / = /i + ji with measurable, nonnegative functions /i,/2 
such that ||/i||i := j fi(x)dx, \\f2W1 > 0. Assume that Y\, I2 are random 
variables with densities /i/||/i||i and /2/H/2II1 respectively and that B is a 
Bernoulli ( 1 1 /1 1 1 1 ) random variable independent of (Yi, Y2). Then the random 
variable BY\ + (1 — -B)Y2 has density /. 

The aim now is to obtain for the densities <p x of Ux + U(1 — U) represen- 
tations ip x = r + g x as above, where r is independent of x G [0, 1]. Typically 
this may not be possible since one may have inf ip x = such that a non-zero 
r independent of x does not exist. However, in our particular situation we 



have Q, hence we are able to choose, e.g., 

r(t) := il [0)1 / 4 ) (t), t€[0,l]. (3) 

Clearly, C//4 has density r/||r||i and let us assume for the moment that a ran- 
dom variable Y x with density 5x/||5a;||i can be simulated via its inverse dis- 
tribution function (quantile function) G~ , i.e., C(Y X ) = C(G~ 1 (U)). Then, 
with a Bernoulli(||r||i) random variable B, independent of U, we have that 

for all x € [0, 1] 

Ux + U(l-U)±^ + (l-B)GZ\U). 

Hence, our update function is $ : [0, 1] x {0, 1} x [0, 1] — > [0, 1], (x, b, u) i-> 
bu/A + (1 — b)G~ l {u). If we construct our Markov chain from the past 
using <I> , in each step there is a probability of ||r||i = 1/8 that all chains 
couple simultaneously. In other words, we can just start at a Geometric (1/8) 
distributed time N in the past, the first instant of {B = 1} when moving 
back into the past. At this time — N we couple all chains via X_^r := U-n/4 
and let the chain run from there until time using the updates GZ-.(Uj+i) 
for j = —N, . . . , —1. It is shown in |10[ Section 2.1] that this is a valid 
implementation of the coupling from the past algorithm in general. 

Hence, we need to derive expressions for the functions G" 1 containing 
only elementary operations. It was calculated in [91 equation (28)] that, for 
all t £ [0, 1] we have 

?*(*) = ((i + *o 2 - 4 <r 1/2 ( w*) + 2 • i [xM (t)) 



with b x := ((1 + x)/2) 2 . Hence, with r given in O) we have (p x (t) > r{t) for 
all x, t & [0, 1]. Note that coupling occurs faster when the function r can be 
chosen larger. For our densities ip x we could as well choose 

r*(t) = — ±= l [om (t), t€[0,l]. 

Then we have <p x (t) > r*(i) > r(t) for all x,t E [0, 1]. However, the subse- 
quent inversion of distribution functions can be done elementary with our 
choice of r. 

We need to invert the distribution functions G x : [0, 1] — > [0, 1] corre- 
sponding to the normalized versions of g x = ip x — r. We have 

gm = V \ Y - r v*® - r ^ dt= ^ i F *^ ~i( yA v< 

!- Inli Jo 7 v 2 



where 



F x {y) ■■- 



\ 1 + x- y/(l + x) 2 



Ay 



l, 



%, 



< y < x, 

x <y <b x , 
b x <y<l, 



is the distribution function of Ux + ^(1 — U). 

The inversion of G x can be done by explicit calculations and yields 

-\z + y/7z+ {I- x) 2 + x -I, ifx€[0,V4],z€[0,g x ], 



G-\z) 
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+ 2 v / 7z + 9 + x(x + 2)-6, ifxe[0,i/4],zGfe,rJ, 



2^(15 + 8x - 7z)(l + 8x + 7z), if x G [0, 1/4], « G (r x , 1], 

-|z + v /7z+(l-x) 2 + x-l, ifxG(V4,l],^G[0, Scc ], 

^(7 + 8x - 7*)(1 + 7z), if x G (1/4, 1], z G (s x , t x ], 

2^(15 + 8x - 7z)(l + 8x + 7z), if x G (1/4, 1], z G (t x , 1], 



where 



l--Vx(x + 2), 



Sx := y (3 + 4x - 4Vx(x + 2)) , ^ := -(8x - 1). 

3 The algorithm 

Our algorithm Simulate [Y =d UY+U(1 — U)~\ has the form discussed in the 
previous section: It draws back to a sequence of independent uniform[0, 1] 
random variables {U- n )n>0 and an independent geometrically distributed 
random variable. (This clearly can be simulated on the basis of independent 
uniform [0, 1] random variables as well.) 

Simulate [Y = d UY + U{\ - [/)] : 



N i- Geometric(l/8) 

X <- U- N /4 

for j from — N + 1 to do X 

return(X) 



GxVj) 



Ilk. 



Figure 1: Histogram of the values of 10 million independent samples from 
/j, generated with the algorithm Simulate [Y =d UY + (7(1 — U)~\ . 

The analysis of the complexity of this algorithm is trivial as the loop is 
iterated a random Geometric(l/8) number of times, hence, e.g., on average 
eight times. 

In Figure 1 the histogram (normalized to area 1) of the values of 10 
million independent samples generated with Simulate \Y =d UY+U(1—U)'] 
is plotted. This simulation was done within a few seconds. A numerical 
approximation of the density of [i has already been presented in [9, Figure 
1]. 
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